function y = func_delta_m0(CONSTS, plot_data)
    
    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eta   = CONSTS.eta;
    eps_a = CONSTS.eps_a;
    
    period = pi/(k0*a);
    dots_per_period = 100;
    num_periods = 3;
    p = (1.0 : period/dots_per_period : num_periods*period)';


    s      = func_s(p, CONSTS);
    n1     = func_n1_from_p(p, CONSTS);
    n2     = func_n2_from_p(p, CONSTS);
    q1     = func_q1_from_p(p, CONSTS);
    q2     = func_q2_from_p(p, CONSTS);

    y = -(n2-n1).*((1./((k0*a*s).^2)) + (eta/eps_a)*(1./(((k0*a)^2)*q1.*q2))) - ...
        (1./(((k0*a)^2)*q1.*s)).*(n2-(eta/eps_a)*n1) - ...
        (1./(((k0*a)^2)*q2.*s)).*((eta/eps_a)*n2-n1);
    
    if (plot_data)
        figure; plot(p, real(y), 'b.-'); title('real(determinant), m=0');
                 xlabel('p'); ylabel('\Delta_0');
        figure; plot(p, imag(y), 'b.-'); title('imag(determinant), m=0');
                 xlabel('p'); ylabel('\Delta_0');
    end

end

function y = func_s(p, CONSTS)
    eps_a = CONSTS.eps_a;
    y = (p.^2 - eps_a).^0.5;
end
